Hypothesis: We hypothesize that at the present time the northern bald ibis population can survive without further management and release. We predict that the observed demographic rates will ensure population growth and do not differ between the breeding colonies.
Study area: Austria, Germany and Italy
Data: Data of the mean per scenario for the first NetLogo model

In this script, we created GAMs and GLMs to find out whether survival per stage or reproduction probabilities have a significant effect on lambda.

Setup

## for non-CRAN packages please keep install instruction
## but commented so it is not run each time, e.g.
# devtools::install_github("EcoDynIZW/template")

## libraries used in this script
## please add ALL LIBRARIES NEEDED HERE
## please remove libraries from the list that are not needed anymore 
## at a later stage
library("effects")

## Loading required package: carData

## lattice theme set by effectsTheme()
## See ?effectsTheme for details.

library("gam")

## Loading required package: splines

## Loading required package: foreach

## Loaded gam 1.20

library("ggeffects")
library("here")

## here() starts at C:/Users/sinah/Documents/IZW/Drenske_2020_PVA_NBI

library("lattice")
library("mgcv")

## Loading required package: nlme

## This is mgcv 1.8-33. For overview type 'help("mgcv-package")'.

## 
## Attaching package: 'mgcv'

## The following objects are masked from 'package:gam':
## 
##     gam, gam.control, gam.fit, s

library("modelsummary")
library('MuMIn')
library("plot3D")

## Warning: package 'plot3D' was built under R version 4.0.5

library("pwr")
library("tidyverse")

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.2
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.4.0     v forcats 0.5.0

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x purrr::accumulate() masks foreach::accumulate()
## x dplyr::collapse()   masks nlme::collapse()
## x dplyr::filter()     masks stats::filter()
## x dplyr::lag()        masks stats::lag()
## x purrr::when()       masks foreach::when()

Data

df_sc_mean <- readRDS(file = here("output", "data-proc", "model-proc", "df_sc_mean.rds"))

Management improvement scenarios

unique(df_sc_mean$Repro)

## [1] 0.53 0.58 0.66 1.06 1.41 3.97

Change the column names

colnames(df_sc_mean)[20:24] <- c("m1", "m2", "m3", "m4", "RR")

Add columns for survival

df_sc_mean$s1 <- (1 - df_sc_mean$m1)
df_sc_mean$s2 <- (1 - df_sc_mean$m2)
df_sc_mean$s3 <- (1 - df_sc_mean$m3)
df_sc_mean$s4 <- (1 - df_sc_mean$m4)

Mortality and respective survival

print("m1 vs s1")

## [1] "m1 vs s1"

table(df_sc_mean$m1)

## 
##  0.2  0.3 0.36 
##  108  108  110

table(df_sc_mean$s1)

## 
## 0.64  0.7  0.8 
##  110  108  108

print("m2 vs s2")

## [1] "m2 vs s2"

table(df_sc_mean$m2)

## 
## 0.08 0.19 0.26 
##  108  108  110

table(df_sc_mean$s2)

## 
## 0.74 0.81 0.92 
##  110  108  108

print("m3 vs s3")

## [1] "m3 vs s3"

table(df_sc_mean$m3)

## 
## 0.14 0.24 0.31 
##  108  108  110

table(df_sc_mean$s3)

## 
## 0.69 0.76 0.86 
##  110  108  108

print("m4 vs s4")

## [1] "m4 vs s4"

table(df_sc_mean$m4)

## 
## 0.02 0.14 0.22 
##  108  108  110

table(df_sc_mean$s4)

## 
## 0.78 0.86 0.98 
##  110  108  108

Calculate extinction probability

df_sc_mean$pext <- df_sc_mean$Extinct/100

Explore database

names(df_sc_mean)

##  [1] "Scenario"        "Run"             "Year"            "IndAll"          "IndB"            "IndK"            "IndUb"           "IndFP"           "IndP"           
## [10] "IndBFP"          "IndBP"           "IndKFP"          "IndKP"           "IndUbFP"         "IndUbP"          "IndJuv"          "IndSub1"         "IndSub2"        
## [19] "IndAdu"          "m1"              "m2"              "m3"              "m4"              "RR"              "N_Supplements"   "Supplement_Time" "Stoch_Frequency"
## [28] "Stoch_Severity"  "N_Juveniles"     "N_Subadults1"    "N_Subadults2"    "N_Adults"        "Lambda_sto"      "SD_Lambda"       "Extinct"         "s1"             
## [37] "s2"              "s3"              "s4"              "pext"

View(df_sc_mean) # okay, these are mortalities.

unique(df_sc_mean[,c(20:24)]) # m1, m2, m3, m4, RR values per scenario

## # A tibble: 326 x 5
##       m1    m2    m3    m4    RR
##    <dbl> <dbl> <dbl> <dbl> <dbl>
##  1   0.2  0.08  0.14  0.02 0.53 
##  2   0.2  0.08  0.14  0.02 0.580
##  3   0.2  0.08  0.14  0.02 0.66 
##  4   0.2  0.08  0.14  0.02 1.06 
##  5   0.2  0.08  0.14  0.14 0.53 
##  6   0.2  0.08  0.14  0.14 0.580
##  7   0.2  0.08  0.14  0.14 0.66 
##  8   0.2  0.08  0.14  0.14 1.06 
##  9   0.2  0.08  0.14  0.22 0.53 
## 10   0.2  0.08  0.14  0.22 0.580
## # ... with 316 more rows

unique(df_sc_mean$s1)

## [1] 0.80 0.70 0.64

unique(df_sc_mean$s2)

## [1] 0.92 0.81 0.74

unique(df_sc_mean$s3)

## [1] 0.86 0.76 0.69

unique(df_sc_mean$s4)

## [1] 0.98 0.86 0.78

unique(df_sc_mean$RR)

## [1] 0.53 0.58 0.66 1.06 1.41 3.97

# start with pExt
table(df_sc_mean$pext)#df_sc_StevS_mean

## 
##    0 0.01 0.02 0.03 0.06 0.24 
##  316    4    3    1    1    1

#  0    0.01 0.02 0.03 0.06 0.24 
# 316    4    3    1    1    1 

# Comment: Only in 10 out of 326 scenarios did the population ever go extinct. 
# With such a distribution of the response variable, 
# there is no need to calculate a model. 
# Therefore we have described the distribution of survival and reproduction values

rows_ext_scen <- which(df_sc_mean$pext > 0)
ext_scen <- df_sc_mean[rows_ext_scen,]
ext_scen

## # A tibble: 10 x 40
##    Scenario    Run  Year IndAll  IndB  IndK IndUb IndFP  IndP IndBFP IndBP IndKFP IndKP IndUbFP IndUbP IndJuv IndSub1 IndSub2 IndAdu    m1    m2    m3    m4    RR N_Supplements
##    <fct>     <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl>   <dbl>  <dbl>  <dbl>   <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>         <dbl>
##  1 0.20.26~ 10451.  24.9   47.4     0     0     0     0  41.3      0     0      0     0       0      0  0.728   11.7     8.67   26.3  0.2   0.26  0.31  0.22 0.53              0
##  2 0.30.19~ 17651.  25.0   42.7     0     0     0     0  36.6      0     0      0     0       0      0  0.727    9.65    7.83   24.5  0.3   0.19  0.31  0.22 0.53              0
##  3 0.30.26~ 21250.  24.9   31.7     0     0     0     0  25.8      0     0      0     0       0      0  0.729    7.37    5.51   18.1  0.3   0.26  0.31  0.22 0.53              0
##  4 0.30.26~ 21351.  25.0   41.4     0     0     0     0  35.5      0     0      0     0       0      0  0.727    9.98    7.43   23.3  0.3   0.26  0.31  0.22 0.580             0
##  5 0.360.1~ 27250.  25     42.9     0     0     0     0  36.7      0     0      0     0       0      0  0.725    9.18    7.46   25.5  0.36  0.19  0.24  0.22 0.53              0
##  6 0.360.1~ 28451.  24.9   31.9     0     0     0     0  26.0      0     0      0     0       0      0  0.730    6.90    5.71   18.6  0.36  0.19  0.31  0.22 0.53              0
##  7 0.360.1~ 28551.  24.9   42.1     0     0     0     0  36.1      0     0      0     0       0      0  0.727    9.52    7.76   24.1  0.36  0.19  0.31  0.22 0.580             0
##  8 0.360.2~ 30851.  24.9   34.2     0     0     0     0  28.2      0     0      0     0       0      0  0.729    7.51    5.69   20.2  0.36  0.26  0.24  0.22 0.53              0
##  9 0.360.2~ 50051.  24.3   21.4     0     0     0     0  15.6      0     0      0     0       0      0  0.750    4.84    3.54   12.3  0.36  0.26  0.31  0.22 0.53              0
## 10 0.360.2~ 32151.  25.0   29.8     0     0     0     0  24.1      0     0      0     0       0      0  0.727    6.90    5.21   17.0  0.36  0.26  0.31  0.22 0.580             0
## # ... with 15 more variables: Supplement_Time <dbl>, Stoch_Frequency <dbl>, Stoch_Severity <dbl>, N_Juveniles <dbl>, N_Subadults1 <dbl>, N_Subadults2 <dbl>, N_Adults <dbl>,
## #   Lambda_sto <dbl>, SD_Lambda <dbl>, Extinct <int>, s1 <dbl>, s2 <dbl>, s3 <dbl>, s4 <dbl>, pext <dbl>

table(ext_scen$m1)

## 
##  0.2  0.3 0.36 
##    1    3    6

table(ext_scen$m2)

## 
## 0.19 0.26 
##    4    6

table(ext_scen$m3)

## 
## 0.24 0.31 
##    2    8

table(ext_scen$m4)

## 
## 0.22 
##   10

table(ext_scen$RR)

## 
## 0.53 0.58 
##    7    3

# We can see that the highest mortality of adults was part of all scenarios where the population went extinct
# See also the excel file under "output\data-proc\model-proc\08_PVA_GLM_table_ext_scen.xlsx"

# Explore lambda for s1 and RR
x <- df_sc_mean$s1
y <- df_sc_mean$RR
z <- df_sc_mean$Lambda_sto
scatter3D(x, y, z, pch = 19, cex = 0.7, clab = c("Lambda"), xlab = "s1", ylab = "RR", zlab = "Lambda")

contourplot(Lambda_sto ~ s1 * RR, data = df_sc_mean, labels = FALSE,
            region = TRUE)

plot(df_sc_mean$Lambda_sto)
abline(h=1) # Lambda > 1 indicates population growth

# Explore lambda for s4 and RR
x <- df_sc_mean$s4
y <- df_sc_mean$RR
z <- df_sc_mean$Lambda_sto
scatter3D(x, y, z, pch = 19, cex = 0.7, clab = c("Lambda"), xlab = "s4", ylab = "RR", zlab = "Lambda")

contourplot(Lambda_sto ~ s4 * RR, data = df_sc_mean, labels = FALSE,
            region = TRUE)

GAM

A general additive model (GAM) tests for non-linearity

The extinction probability is the response variable and the predictor variables are formed by the survival probabilities per stage (s1-s4) and the reproduction (RR)

# define a k
# myk - stands for 'my number of knots/ smoothing splines' :-)
myk <- 3 # or 4

Response = Lambda

## Lambda and survival
gammod_2_1 <- gam(formula = Lambda_sto ~ s(s1, k=myk) + s(s2, k=myk) + 
                          s(s3, k=myk) + s(s4, k=myk) +  
                          s(RR, k=myk), 
                          data = df_sc_mean, 
                          na.action = 'na.fail', family=Gamma) 
summary(gammod_2_1)

## 
## Family: Gamma 
## Link function: inverse 
## 
## Formula:
## Lambda_sto ~ s(s1, k = myk) + s(s2, k = myk) + s(s3, k = myk) + 
##     s(s4, k = myk) + s(RR, k = myk)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.9036404  0.0004053    2229   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df      F p-value    
## s(s1) 1.000  1.000  964.7  <2e-16 ***
## s(s2) 1.000  1.000  887.0  <2e-16 ***
## s(s3) 1.000  1.000  893.8  <2e-16 ***
## s(s4) 1.908  1.991 5617.9  <2e-16 ***
## s(RR) 1.999  2.000 5328.9  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.985   Deviance explained = 98.6%
## GCV = 6.6955e-05  Scale est. = 6.5312e-05  n = 326

# F values for s4 and RR are very high - they have the highest influence

## Lambda and mortality
gammod_2_2 <- gam(formula = Lambda_sto ~ s(m1, k=myk) + s(m2, k=myk) + 
                          s(m3, k=myk) + s(m4, k=myk) +  
                          s(RR, k=myk), 
                          data = df_sc_mean, 
                          na.action = 'na.fail', family=Gamma) 
summary(gammod_2_2)

## 
## Family: Gamma 
## Link function: inverse 
## 
## Formula:
## Lambda_sto ~ s(m1, k = myk) + s(m2, k = myk) + s(m3, k = myk) + 
##     s(m4, k = myk) + s(RR, k = myk)
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.9036404  0.0004053    2229   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##         edf Ref.df      F p-value    
## s(m1) 1.000  1.000  964.7  <2e-16 ***
## s(m2) 1.000  1.000  887.0  <2e-16 ***
## s(m3) 1.000  1.000  893.8  <2e-16 ***
## s(m4) 1.908  1.991 5617.9  <2e-16 ***
## s(RR) 1.999  2.000 5328.9  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.985   Deviance explained = 98.6%
## GCV = 6.6955e-05  Scale est. = 6.5312e-05  n = 326

# F values for s4 and RR are very high - they have the highest influence

# no difference between survival and mortality

Plots

plot(gammod_2_1)

vis.gam(gammod_2_1, view = c("s4", "RR"), plot.type = "persp", type = "response", zlim = c(0.5,1.8), 
        ticktype = "detailed", border = NA, n.grid = 1000, nticks = 3, 
        xlab = "s4", ylab = "RR", zlab = "Lambda", 
        color = "heat", theta= -140, cex.lab = 3, cex.axis = 2.5)

vis.gam(gammod_2_2, view = c("m4", "RR"), plot.type = "persp", type = "response", zlim = c(0.5,1.8), 
        ticktype = "detailed", border = NA, n.grid = 1000, nticks = 3, 
        xlab = "m4", ylab = "RR", zlab = "Lambda", 
        color = "heat", theta= -140, cex.lab = 3, cex.axis = 2.5)

Save the plot

# # s4 and RR
# png(filename = here("plots", "04_PVA", "persp_s4_RR_lambda_20211124.png"), width = 1200, height = 772)
# vis.gam(gammod_2_1, view = c("s4", "RR"), plot.type = "persp", type = "response", zlim = c(0.5,1.8), 
#         ticktype = "detailed", border = NA, n.grid = 1000, nticks = 3, 
#         xlab = "s4", ylab = "RR", zlab = "Lambda", 
#         color = "heat", theta= -140, cex.lab = 3, cex.axis = 2.5)
# dev.off()

Effect size

Calculating the effect size

pwr.t.test(n=326, sig.level=0.05, power = 0.8, alternative="greater")

## 
##      Two-sample t test power calculation 
## 
##               n = 326
##               d = 0.1949595
##       sig.level = 0.05
##           power = 0.8
##     alternative = greater
## 
## NOTE: n is number in *each* group

GLM

Response = Lambda

# With interactions between RR and all survival values
mod1_3_1 <- glm(Lambda_sto ~ (s1 + s2 + s3 + s4) + RR, 
              data = df_sc_mean, na.action = 'na.fail', family=Gamma)

# With interactions between RR and all survival values
mod1_3_2 <- glm(Lambda_sto ~ (m1 + m2 + m3 + m4) + RR, 
              data = df_sc_mean, na.action = 'na.fail', family=Gamma)

summary(mod1_3_1)

## 
## Call:
## glm(formula = Lambda_sto ~ (s1 + s2 + s3 + s4) + RR, family = Gamma, 
##     data = df_sc_mean, na.action = "na.fail")
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.174754  -0.010237  -0.002538   0.005698   0.039817  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.870086   0.018618  100.45   <2e-16 ***
## s1          -0.199206   0.012780  -15.59   <2e-16 ***
## s2          -0.170727   0.011390  -14.99   <2e-16 ***
## s3          -0.182280   0.012103  -15.06   <2e-16 ***
## s4          -0.522307   0.010212  -51.15   <2e-16 ***
## RR          -0.122401   0.002476  -49.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.000289145)
## 
##     Null deviance: 1.507876  on 325  degrees of freedom
## Residual deviance: 0.095412  on 320  degrees of freedom
## AIC: -1645.8
## 
## Number of Fisher Scoring iterations: 4

anova(mod1_3_1, test='LRT')

## Analysis of Deviance Table
## 
## Model: Gamma, link: inverse
## 
## Response: Lambda_sto
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                   325    1.50788              
## s1    1  0.05411       324    1.45377 < 2.2e-16 ***
## s2    1  0.04857       323    1.40520 < 2.2e-16 ***
## s3    1  0.04808       322    1.35712 < 2.2e-16 ***
## s4    1  0.67890       321    0.67822 < 2.2e-16 ***
## RR    1  0.58281       320    0.09541 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# s4 and RR have the highest deviance - they explain the model the most

#dredge(mod1_3_1)
plot(allEffects(mod1_3_1),type='response', ylim = c(1,2))

Save the plot

png(filename = here("plots", "04_PVA", "lambda_effects_plot20211216.png"), width = 1200, height = 772)
plot(allEffects(mod1_3_1),type='response', ylim = c(1,2))
dev.off()

## png 
##   2

R square R^2= 1-(Residual deviance/Null deviance)

round(1 - (0.095412/1.507876), digits = 2)

## [1] 0.94

summary(mod1_3_2)

## 
## Call:
## glm(formula = Lambda_sto ~ (m1 + m2 + m3 + m4) + RR, family = Gamma, 
##     data = df_sc_mean, na.action = "na.fail")
## 
## Deviance Residuals: 
##       Min         1Q     Median         3Q        Max  
## -0.174754  -0.010237  -0.002538   0.005698   0.039817  
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.795567   0.005226  152.23   <2e-16 ***
## m1           0.199206   0.012780   15.59   <2e-16 ***
## m2           0.170727   0.011390   14.99   <2e-16 ***
## m3           0.182280   0.012103   15.06   <2e-16 ***
## m4           0.522307   0.010212   51.15   <2e-16 ***
## RR          -0.122401   0.002476  -49.43   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.000289145)
## 
##     Null deviance: 1.507876  on 325  degrees of freedom
## Residual deviance: 0.095412  on 320  degrees of freedom
## AIC: -1645.8
## 
## Number of Fisher Scoring iterations: 4

anova(mod1_3_2, test='LRT')

## Analysis of Deviance Table
## 
## Model: Gamma, link: inverse
## 
## Response: Lambda_sto
## 
## Terms added sequentially (first to last)
## 
## 
##      Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                   325    1.50788              
## m1    1  0.05411       324    1.45377 < 2.2e-16 ***
## m2    1  0.04857       323    1.40520 < 2.2e-16 ***
## m3    1  0.04808       322    1.35712 < 2.2e-16 ***
## m4    1  0.67890       321    0.67822 < 2.2e-16 ***
## RR    1  0.58281       320    0.09541 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# s4 and RR have the highest deviance - they explain the model the most

#dredge(mod1_3_2)
plot(allEffects(mod1_3_2),type='response')

R square R^2= 1-(Residual deviance/Null deviance)

round(1 - (0.095412/1.507876), digits = 2)

## [1] 0.94

Model summary

datasummary_skim(df_sc_mean)

	Unique (#)	Mean	SD	Min	Median	Max
Run	326	16525.0	10085.5	50.5	16300.5	60150.5
Year	207	19.5	5.6	6.6	20.3	25.0
IndAll	326	830.9	473.9	21.4	1162.3	1406.3
IndB	1	0.0	0.0	0.0	0.0	0.0
IndK	1	0.0	0.0	0.0	0.0	0.0
IndUb	1	0.0	0.0	0.0	0.0	0.0
IndFP	1	0.0	0.0	0.0	0.0	0.0
IndP	326	813.2	465.2	15.6	1143.5	1386.3
IndBFP	1	0.0	0.0	0.0	0.0	0.0
IndBP	1	0.0	0.0	0.0	0.0	0.0
IndKFP	1	0.0	0.0	0.0	0.0	0.0
IndKP	1	0.0	0.0	0.0	0.0	0.0
IndUbFP	1	0.0	0.0	0.0	0.0	0.0
IndUbP	1	0.0	0.0	0.0	0.0	0.0
IndJuv	193	1.0	0.4	0.7	0.9	2.6
IndSub1	326	214.3	135.0	4.8	252.5	674.8
IndSub2	326	152.1	91.8	3.5	181.6	351.3
IndAdu	326	463.5	266.3	12.3	576.7	787.2
m1	3	0.3	0.1	0.2	0.3	0.4
m2	3	0.2	0.1	0.1	0.2	0.3
m3	3	0.2	0.1	0.1	0.2	0.3
m4	3	0.1	0.1	0.0	0.1	0.2
RR	6	0.7	0.3	0.5	0.7	4.0
N_Supplements	1	0.0	0.0	0.0	0.0	0.0
Supplement_Time	1	4.0	0.0	4.0	4.0	4.0
Stoch_Frequency	1	0.0	0.0	0.0	0.0	0.0
Stoch_Severity	1	0.2	0.0	0.2	0.2	0.2
N_Juveniles	1	37.0	0.0	37.0	37.0	37.0
N_Subadults1	1	11.0	0.0	11.0	11.0	11.0
N_Subadults2	1	8.0	0.0	8.0	8.0	8.0
N_Adults	1	18.0	0.0	18.0	18.0	18.0
Lambda_sto	326	1.1	0.1	1.0	1.1	1.4
SD_Lambda	326	0.0	0.0	0.0	0.0	0.0
Extinct	6	0.1	1.4	0	0.0	24
s1	3	0.7	0.1	0.6	0.7	0.8
s2	3	0.8	0.1	0.7	0.8	0.9
s3	3	0.8	0.1	0.7	0.8	0.9
s4	3	0.9	0.1	0.8	0.9	1.0
pext	6	0.0	0.0	0.0	0.0	0.2

varnam <- c('s1' = 'Survival of stage 1',
            's2' = 'Survival of stage 2',
            's3' = 'Survival of stage 3',
            's4' = 'Survival of stage 4',
            'RR' = 'Reproductive Rate')
varnam2 <- c('m1' = 'Mortality of stage 1',
            'm2' = 'Mortality of stage 2',
            'm3' = 'Mortality of stage 3',
            'm4' = 'Mortality of stage 4',
            'RR' = 'Reproductive Rate')


# modelplot(mod1_3,
#           coef_omit = 'Interc',
#           coef_map = varnam,
#           conf_level = .99) 

modelplot(mod1_3_1, 
          coef_omit = 'Interc',
          coef_map = varnam, 
          conf_level = .99)

modelplot(mod1_3_2, 
          coef_omit = 'Interc',
          coef_map = varnam2, 
          conf_level = .99)

# You can see that s4/m4 and RR are important and s1-s3/m1-m3 are not.

Session Info

## DO NOT REMOVE!
## We store the settings of your computer and the current versions of the
## packages used to allow for reproducibility
Sys.time()

## [1] "2022-01-06 18:55:07 CET"

#git2r::repository() ## uncomment if you are using GitHub
sessionInfo()

## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                    LC_TIME=German_Germany.1252    
## 
## attached base packages:
## [1] splines   stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.0      stringr_1.4.0      dplyr_1.0.2        purrr_0.3.4        readr_1.4.0        tidyr_1.1.2        tibble_3.0.3       ggplot2_3.3.2      tidyverse_1.3.0   
## [10] pwr_1.3-0          plot3D_1.4         MuMIn_1.43.17      modelsummary_0.6.5 mgcv_1.8-33        nlme_3.1-149       lattice_0.20-41    here_0.1           ggeffects_1.0.0   
## [19] gam_1.20           foreach_1.5.1      effects_4.2-0      carData_3.0-4     
## 
## loaded via a namespace (and not attached):
##  [1] minqa_1.2.4       colorspace_1.4-1  ellipsis_0.3.1    showtext_0.9      sjlabelled_1.1.7  rprojroot_1.3-2   estimability_1.3  fs_1.5.0          rstudioapi_0.11  
## [10] farver_2.0.3      showtextdb_3.0    remotes_2.2.0     fansi_0.4.1       lubridate_1.7.9   xml2_1.3.2        codetools_0.2-16  knitr_1.30        pkgload_1.1.0    
## [19] jsonlite_1.7.1    nloptr_1.2.2.2    broom_0.7.4       dbplyr_1.4.4      compiler_4.0.3    d6_0.1.0.0        httr_1.4.2        backports_1.1.10  assertthat_0.2.1 
## [28] Matrix_1.2-18     survey_4.0        cli_2.0.2         htmltools_0.5.0   prettyunits_1.1.1 tools_4.0.3       misc3d_0.9-0      gtable_0.3.0      glue_1.4.2       
## [37] tables_0.9.6      tinytex_0.26      Rcpp_1.0.5        cellranger_1.1.0  vctrs_0.3.4       iterators_1.0.13  insight_0.11.1    xfun_0.18         ps_1.4.0         
## [46] testthat_2.3.2    lme4_1.1-25       rvest_0.3.6       lifecycle_0.2.0   devtools_2.3.2    statmod_1.4.35    MASS_7.3-53       scales_1.1.1      hms_0.5.3        
## [55] yaml_2.2.1        memoise_1.1.0     stringi_1.5.3     highr_0.8         desc_1.2.0        checkmate_2.0.0   boot_1.3-25       pkgbuild_1.1.0    rlang_0.4.7      
## [64] pkgconfig_2.0.3   evaluate_0.14     labeling_0.3      processx_3.4.4    tidyselect_1.1.0  magrittr_1.5      bookdown_0.20     R6_2.4.1          generics_0.0.2   
## [73] DBI_1.1.0         pillar_1.4.6      haven_2.3.1       withr_2.3.0       survival_3.2-7    nnet_7.3-14       modelr_0.1.8      crayon_1.3.4      utf8_1.1.4       
## [82] rmarkdown_2.4     sysfonts_0.8.1    usethis_1.6.3     grid_4.0.3        readxl_1.3.1      blob_1.2.1        callr_3.5.0       rmdformats_0.3.7  webshot_0.5.2    
## [91] reprex_0.3.0      digest_0.6.25     stats4_4.0.3      munsell_0.5.0     viridisLite_0.3.0 kableExtra_1.3.1  mitools_2.4       tcltk_4.0.3       sessioninfo_1.1.1

PVA of the NBI: GLM of survival and reproduction probabilities

PVA of the NBI: GLM of survival and reproduction probabilities

Setup

Data

Management improvement scenarios

Calculate extinction probability

Explore database

GAM

Response = Lambda

Plots

Effect size

GLM

Response = Lambda

Model summary